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ABSTRACT 

We utilise ground-based, balloon-borne and satellite climatology data to reconstruct 
site and season-dependent vertical profiles of precipitable water vapour (PWV). We 
use these profiles to solve radiative transfer through the atmosphere, and derive at¬ 
mospheric brightness temperature (Tatm) and optical depth (t) at centimetre wave¬ 
lengths. 

We validate the reconstruction by comparing the model column PWV with photo¬ 
metric measurements of PWV, performed in clear sky conditions pointed towards the 
Sun. Based on the measurements, we devise a selection criteria to filter the climatology 
data to match the PWV levels to the expectations of the clear sky conditions. 

We apply the reconstruction to the location of a Polish 32-metre radio telescope, 
and characterise Tatm and r year-round, at selected frequencies. We also derive the 
zenith distance dependence for these parameters, and discuss the shortcomings of using 
planar, single-layer, and optically thin atmospheric models in continuum radio-source 
flux-density measurement calibrations. 

We obtain PWV-Tatm and PWV-r scaling relations in clear sky conditions, and 
constrain limits to which the actual Tatm and r can deviate from those derived solely 
from the climatological data. 

Finally, we suggest a statistical method to detect clear sky that involves ground- 
level measurements of relative humidity. Accuracy is tested using local climatological 
data. The method may be useful to constrain cloud cover in cases when no other (and 
more robust) climatological data are available. 

Key words: radiative transfer - atmospheric effects - site testing - radio continuum: 
general - methods: observational 


1 INTRODUCTION 

Ground-based, cm-wavelengths radio continuum observa¬ 
tions of astrophysical sources depend on atmospheric emis¬ 
sion, and absorption of the incident radiation. Time varying 
line-of-sight (LOS) abundance of ice, liquid water, and wa¬ 
ter vapour generates variable optical depth, which leads to 
signal instabilities. These instabilities do not average out 
under long integrations due to the steep spectrum of tur¬ 
bulent atmospheric water density fluctuations. The atmo¬ 
spheric effects associated with dry air thermal emission, con¬ 
tinuum and line absorption are frequency-dependent, and 
can be characterised by atmospheric brightness temperature 
(Tatm(i^)) and transmittance (tatm(i')). These atmospheric 
emissions contribute to the system temperature (Tsys) that 
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limits the sensitivity of any ground-based telescope-receiver 
pair. 

Over the last few decades, monitoring temperature, 
pressure, density, and other altitude dependent parameters 
of air, have helped to develop a few widely-accepted models 
of Earth’s atmosphere. Advances in atomic and molecular 
line spectroscopy provided absorption coefficients for vari¬ 
ous gas species, and computer-generated spectra for mix¬ 
tures of gases under given thermodynamic conditions can 
now reproduce these observations in great detail. Thus, the 
radiative properties of the atmosphere are well known, and 
can be derived from the first principles, from radio waves to 
infrared frequencies. 

At cm-wavelengths the electromagnetic spectrum of 
Tatm is predominantly defined by continuum absorption, and 
emission of oxygen (O2) and water vapour (H2O) molecules. 
The latter results from strong coupling of the electric dipole 
moment of water molecules to the millimetre radiation via 
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rotational transitions. Going to liquid and ice states less and 
no rotational freedom is possible respectively, and therefore 
droplets and ice particles are expected to radiate less per 
molecule. In the presence of clouds, the background signals 
are additionally attenuated in ice and droplets, decreasing 
the signal-to-noise ratio (SNR) of any radio source. This 
translates to an increased observational time required to de¬ 
tect the same source at the same significance level, as com¬ 
pared to the situation without clouds. Given that the cloud 
cover is non-uniform and variable, the amount of the attenu¬ 
ation will change over time, generating variance at different 
time scales in the signal received through the main beam, 
or through side-lobes. For the radio continuum flux density 
observations at cm- and mm-wavelengths, this means that 
nearly clear sky conditions are required, although in prac¬ 
tice, observations are also viable when some high-level icy 
clouds are present. 

Unlike the 'dry’ component of air, the distribution of at¬ 
mospheric water on Earth is strongly dependent on location. 
For this reason many dedicated radio telescopes are built 
at high altitudes and in dry and/or cold climates (such as 
the Atacama desert or the South Pole). This reduces atmo¬ 
spheric emissions, maximises transmittance, improves ther¬ 
mal stability, and rules out cloud attenuation. At other lo¬ 
cations, atmospheric water variations must be monitored to 
optimally match the observational programme to the current 
weather conditions. In clear sky conditions, measurements of 
precipitable water vapour (PWV) can be useful in working 
out optical depths and estimating atmospheric absorption 
corrections of astronomical radio source flux density mea¬ 
surements. 

Atmospheric water vapour can be measured in a num¬ 
ber of ways, including (i) radio sounding, (ii) atmospheric 
delays of GPS-satellite signals, (iii) sun photometry of water 
lines in near-infrared light, (iv) direct radiometric measure¬ 
ments in water bands (e.g. Liljegren et al. (2001)). In the 
larger time scales, PWV can be modelled statistically using 
climatological data. 

Currently, the publicly-available data from ground- 
based meteorological stations, radio sounding, and satellite 
observations allow reconstruction of the vertical structure 
of atmosphere at particular locations of interest. For ex¬ 
ample, the Integrated Global Radiosonde Archive (IGRA) 
data, which we use in this analysis, provides substantial 
aerial coverage and density worldwide, but only a few 
atmospheric parameters are measured at relatively large 
time intervals (only twice a day). Likewise, many satel¬ 
lite data obtained with the solar occultation method also 
have nearly global coverage, but their sensitivity is largely 
limited in the troposphere due to clouds. However, these 
meteorological data can be used to model the total av¬ 
erage PWV on a month-by-month basis. Whether by di¬ 
rect radiometric observations or by solving radiative transfer 
equations, many radio astronomical observation sites have 
Tatm and fatm Calculated at the desired frequencies and the 
time of year (Ajello et al. 1995; Radford & Holdaway 1998; 
Bussmann et al. 2005; Bustos et al. 2014). In this paper we 
will characterise the atmosphere in these terms for the first 
time for the location of the 32-metre radio telescope (RT32) 
in Poland, and show that a similar approach can be easily 
adapted to virtually any other location. The RT32 is one of 
the European VLBI Network stations, currently capable of 


observations in the L, G, K and Ka frequency bands (see 
Lew et al. (2015) for more RT32 specifications). 

The motivation for investigating local atmosphere 
arises from the planned radio source survey with the 30- 
GHz OCRA-f focal plane, 4-pair, beam-switched receiver 
(Lew et al. 2015). In particular, we seek to improve the at¬ 
mospheric model, previously used for the continuum flux 
density measurement calibrations (Gawronski et al. 2010; 
Peel 2010; Lancaster et al. 2011), which requires estimates 
of the atmospheric optical depth (r) at the source zenith 
distance (za). These estimates are typically obtained by per¬ 
forming tipping scans, with the assumption that Tatm(zd) ~ 
T'atm(O) sec(zd). We investigate the limits of validity of this, 
and other approximations, using radiative transfer in the 
atmosphere in clear sky conditions. 

Radio sounding and satellite data are recorded indepen¬ 
dently of weather conditions, and modelling PWV in clear 
skies requires employing selection criteria that would assure 
that the filtered subset of the data corresponds to the clear 
sky conditions at a particular location and observation time. 
A search for such selection criteria, given the aforementioned 
sparsity and deficiencies of the sounding data, is also one of 
the aims of this paper. 

The structure of this paper is as follows. In section 2 
we discuss the publicly available climatological data used in 
our reconstructions of atmospheric profiles, describe the at¬ 
mospheric model parametrisation and introduce our PWV 
measurements. In section 3 we describe the radiative trans¬ 
fer approach to calculating atmospheric transmittances and 
brightness temperatures. The main results are gathered in 
Sec. 4. Discussion and conclusions are in sections 5 and 6 
respectively. 


2 DATA 

We consider contributions to Tatm and fatm from dry air 
(with standard nitrogen to oxygen proportions), ozone (O 3 ) 
and water vapour (H 2 O), and ignore atmospheric trace 
gases. In the following subsections, we describe the climato¬ 
logical data, used to reconstruct a localised, vertical struc¬ 
ture of the atmosphere, which is a starting point for radiative 
transfer. 


2.1 Vertical pressure and temperature profiles 

We reconstruct vertical atmospheric temperature and pres¬ 
sure profiles using the GOSPAR^ International Reference 
Atmosphere project data^ (hereafter GIRA 86 ) (Rees 1988; 
Rees et al. 1990; Rees 1992). These zonally-averaged data 
are a compilation of ground-based, radiosonde and satel¬ 
lite measurements of atmospheric pressure, temperature, 
wind velocity and geopotential height up to an altitude 
of z = 120 km, with 5° latitude resolution, roughly 2 km 
vertical resolution, and nearly global coverage in latitudes 
{9 = [80 °S, 80 °N]). We use the monthly averaged vertical 
profiles of pressure P{zg) and temperature T{zg), given as a 
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Figure 1. Atmospheric pressure profiles in January (solid, blue) 
and in July (dashed, red) interpolated from the CIRA86 data 
for the latitude Ot = 53.1°N. For comparison, the theoretical 
profile (Eq. 1) is plotted for the assumed values of lapse-rate (Lr), 
pressure (Pq), and temperature (To) (bottom solid-green curve). 
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The interpolated CIRA86 pressure profiles are shown 
in Fig. 1. For visualisation we compare them against the 
theoretical model with the altitude-independent lapse rate 
(NASA 1976): 


P{z)=Po[ 


To 




To + Lr{z — Zo} 


( 1 ) 


where Lr is the lapse rate (Fig. 1), where go ~ 
9.80665 km/s^ is the assumed surface acceleration of the 
Earth, R = 8.31432 J ■ mol“^ • is the gas constant, 
g = 0.0289644 kg ■ mol~^ is the molar mass of the air and 
Po, To and zo are fiducial pressure, temperature and altitude 
respectively. While the consistency is very good in the tro¬ 
posphere, the deviation from the measured profiles in the 
stratosphere results from the lapse rate increasing to zero 
as one approaches tropopause, the fact that we do not take 
into account in this visualisation. 

In order to adjust the pressure profiles to the local 
conditions and improve the reconstruction at low altitudes 
{z < 13 km), we combine CIRA86 profile with the profile 
extracted from the local sounding data (see section 2.2.1). 
However, the two profiles turn out to be quite compatible, 
with the biggest discrepancy of 8% at z = 30 km in Jan¬ 
uary, and much smaller in July. 

Similarly, in the case of temperature profiles, we gauge 
the impact of the longitudinal deviations from the CIRA86 
mean, by comparing it with the daytime radio sounding pro¬ 
file from Legionowo (near Warsaw, Poland), averaged over 
the period 2000-2014. We find the two datasets consistent 
to within 8% in the stratosphere and in the lower tropo¬ 
sphere (Fig. 2). The consistency is better between the two 
regimes. The near-ground departures from the zonal aver¬ 
ages are somewhat expected, which is why we rely on the 
sounding data at low altitudes. The two datasets clearly 
mark the transition to tropopause at an altitude of about 
10 km at the considered latitude. 


Figure 2. Atmospheric temperature profiles in January (solid, 
blue) and in July (dashed, red) interpolated from the CIRA86 
data for the latitude Oj- = 53.1°N. The curves extending up to 
z Ri 30 km are the IGRA radiosonde data from Legionowo station, 
averaged over the period of 15 years. The horizontal error bars 
represent Tier dispersion resulting from variable geopotential alti¬ 
tudes assigned to the fixed pressure levels. The vertical error bars 
represent Tier dispersion of temperatures in the selected month. 

function of geopotential height and we apply an interpola¬ 
tion to obtain a profile for the desired latitude. Throughout 
this paper we will interpolate the profiles for the latitude 
9t = 53.1°N - the latitude of the Torun Centre for Astron¬ 
omy (TCfA). Geopotential height (zg) is related to geomet¬ 
rical height (z) above the sea level via: Zg = z R^/{R^ +z), 
where = 6371 km is the assumed value of Earth ra¬ 
dius. In what follows, we will ignore the difference between 
the two heights, as they are almost identical over the range 
of altitudes of interest. The differences are also unimpor¬ 
tant when compared the to variance generated by weather 
related pressure variations. Associating the fixed CIRA86 
pressure levels to altitudes (z) allows us to assign geomet¬ 
rical heights to other physical quantities of interest as well, 
such as temperature, relative humidity, and mixing ratios, 
which are typically recorded at fixed pressure levels. 


2.2 Vertical WV profiles 

2.2.1 Relative Humidity profiles 

Of all trace atmospheric gases, WV is of the greatest rele¬ 
vance for the cm-wavelength observations. WV volume mix¬ 
ing ratio (VMR) — the number density of H 2 O molecules 
relative to number density of all other species in the atmo¬ 
sphere — is highest in the troposphere. At these altitudes 
satellite data are no longer useful, and WV mixing ratios 
or alternatively the relative humidity (RH) for a given pres¬ 
sure and temperature, must be estimated from radio sound¬ 
ing data. The data is publicly available from the Integrated 
Global Radiosonde Archive ® (Durre et al. 2006) updated on 
daily basis. There are several hundreds of stations around 
the globe, allowing vertical profiles reconstruction with the 
land surface spatial resolution of the order a few hundred 
kilometres, on average. A typical radiosonde accuracy to 
perform the temperature, pressure and RH measurement is 
about 0.5 °C, 1 hPa, and 1% respectively (Peixoto & Oort 
1996). 

® http://wwwl.nede.noaa.gov/pub/data/igra 
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Figure 3. {Left) Monthly-average relative humidity profiles from the Legionowo IGRA station in January (triangles) and July (squares) 
for the case of unscreened data (solid), and for the subset of 10% driest conditions (dashed) as indicated by pn = 0.1 in the plot legend. 
{Right) Corresponding H 2 O volume mixing ratio profiles. The vertical line indicates transition altitude above which the satellite data is 
used. The satellite water abundances are not subject to the selection criterion. 


For the purpose of our analysis we choose the ra¬ 
diosonde data from the Legionowo station located near War¬ 
saw (Poland) as they have been archived since 1957. The 
data contain daily records of geopotential height {zg), tem¬ 
perature (T), dew point depression (AT^ = T — Td, where Td 
is dew point temperature) wind direction, and wind speed, 
probed twice a day (midnight and noon) at fixed pressure 
levels {Pi), as the radiosonde travels up the atmosphere. 
We will hereafter refer to this data as IGRA. For any given 
month we use the daily daytime records between and includ¬ 
ing the years 2000 and 2014. We reject the data that have 
incomplete or incorrect height, temperature or dew point 
depression information. This leaves a few thousand records 
per month for further analysis. 

At each pressure level and for each {T,P,Td} tuple we 
derive the corresponding relative humidity by analytically 
solving the following equation: 


Td{T,RR) 


\r{T, RH) 

/? - F(T, RH) 


( 2 ) 


where 


and T is in Celsius and RH is expressed as percentages (see 
e.g. Lawrence (2005) for the derivation). In Eq. 3 the con¬ 
stants are: /? = 17.62 and A = 243.12 °C (Sonntag 1990). 
These values yield Td consistent with a more rigorous deriva¬ 
tion of Hardy (1998) to within 0.3 °C for the temperature 
range from —50 °C to 50 °C, and within the measured range 
of relative humidities. See also Alduchov & Eskridge (1997) 
to review the values of the constants from other studies. 

We associate the RHs with altitudes using the recon¬ 
structed P{z) relation (Sec. 2.1), and employ Akima inter¬ 
polation (Akima 1970) to create a tabulated version of the 
profiles. We use the sounding data up to z 13 km and 
satellite data at higher altitudes (Fig. 3). The near-ground 
tail of each RH profile is supplemented with the average RH 
data point recorded by the local meteorological station. This 
improves adjustment to the RT32 site conditions. 

In order to create a low-humidity subset of the data - 


conditions that typically correspond to a cloudless sky - we 
select the IGRA records according to the quantile function 
Q{Ph) of the RH record distribution, in each of the pres¬ 
sure levels independently (sec. 2.5). Thus pn is a selection 
parameter, corresponding to the probability that a random 
observation at any pressure level will have the RH value 
smaller than Q{pii). In order to visualise the selection effect 
we arbitrarily choose pn = 0.1, corresponding to the 10% 
driest conditions for the considered month in the year-to- 
year data samples (Fig. 3, dashed lines). We will later fine 
tune this choice using external data. The resulting yearly 
samples are averaged into a single profile that approximates 
the local climate for a given month. 

We convert the reconstructed RH profiles into H 2 O 
VMR profiles (a;v(H20)) using the following formula: 


iH nl /RH\Psat(T') 


(4) 


where P is the total atmospheric pressure, and Psat is the 
saturation pressure of water vapour at temperature T. Psat 
is calculated using equation number 10 of Murphy & Koop 
(2005), reported to be valid within the range of temperatures 
considered in this paper. The reconstructed a:v(H20) profiles 
are shown in Fig. 3 (right panel). At higher altitudes the 
profiles are reconstructed using satellite data. 

A comparison of GIRA86 and Legionowo radiosonde 
data in terms of temperature shows that the two agree very 
well in the region of the upper troposphere (Fig. 2). How¬ 
ever, there is a systematical dependence of the upper tro¬ 
pospheric RHs recorded by radiosondes in the period from 
1979 to 1991, on geographical location. The dependence re¬ 
sults from the type of instrumentation that has been used 
(Soden & Lanzante 1996). For this reason, in this analysis, 
we refrain from using data from the previous century, and 
rely on the data from the previous and current decades, as 
technology exchange it is likely to have mitigated these in¬ 
consistencies. Peixoto & Oort (1996) has also investigated 
radiosonde humidity data (recorded in the period from 1973 
to 1988), performed cross-checks with SAGE satellite mea¬ 
surements, and found similar discrepancies as those pointed 
out by Soden & Lanzante (1996), but confirmed that in the 
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z [km] 

Figure 4. WV volume mixing ratio profiles (solid) interpolated 
at 9t = 53.1 °N from ACE-FTS data in January (cyan,thin) and 
July (magenta,thick). The dashed (dash-dotted) lines represent 
the ozone VMR profiles interpolated from the KNMI (ACE-FTS) 
data in January (blue) and July (red) (see Sec. 2.3). 

lower and middle troposphere the discrepancies are not im¬ 
portant. 

2.2.2 Water Vapour mixing ratio profile 

We use the data from the Fourier Transform Spectrometer 
(Jones et al. 2012) of the Atmospheric Chemistry Experi¬ 
ment^ (ACE-FTS) (Bernath et al. 2005) in order to include 
the stratospheric water vapour. 

The interpolated profile is shown in Fig. 4. We combine 
it with the tropospheric radiosonde profile using a linear 
transformation from one profile to another within the range 
of overlapping altitudes (typically 10 < 2 < 13 km). We 
verified, however, that the impact of the stratospheric WV 
is negligible and could be ignored in the case of ground-level 
sites. 

2.3 Ozone 

We model the local atmospheric ozone VMR profile using a 
combination of ozonesonde and satellite data® provided by 
the Royal Netherlands Meteorological Institute (hereafter 
KNMI) The KNMI ozone data are zonally averaged VMRs 
measured over 17 different latitudes from —80°S to 80 °N, 
and at 19 different pressure levels from 0.3 hPa to 1000 hPa. 
The data is compiled from 30 ozone stations and SBUV- 
SBUV/2 satellite observations, collected during the years 
1980-1991(Paul et al. 1998). 

Currently, there is a wealth of ozone data from many 
satellites and from radiosondes (e.g. WOUDC (1961); 
Fioletov et al. (2002); Jones et al. (2012)) that continuously 

^ The Atmospheric Chemistry Experiment (ACE), also known 
as SCISAT, is a Canadian-led mission mainly supported by the 
Canadian Space Agency. 

® http://www.knml.nl/research/climate_chemistry/Data/ 
FKClimatology/ 


monitor the ozone layer. In the current work, our main focus 
is on the cm-wavelengths, where the ozone impact is sub¬ 
dominant, but for the sake of completeness and accounting 
for possible extensions into higher frequencies, we include 
one of the available ozone datasets (KNMI) into the analysis. 
The consistency between the profiles obtained form various 
ozone experiments is shown in Fig. 4. 

2.4 Ground-level meteorological data 

We improve the adjustment of the atmospheric model to 
the local conditions at the ground-level by including data 
recorded by the IRDAM WST7000 meteorological station, 
installed on the roof of TCfA, near RT32. The data include 
temperatures, pressures, and relative humidities recorded at 
a high time resolution. We used the data covering a pe¬ 
riod that slightly exceeds 4 years i.e. from October 2010 to 
November 2014. In each month they contribute as a single 
data point in the RH profile (Fig. 3) at the altitude of TCfA, 
but also influence higher altitudes via smooth interpolation. 
We detect a systematical effect in this RH dataset, resulting 
from a progressive deterioration of the humidity sensor over 
time. The effect leads to an overestimation of the local RHs 
by < 4% in the years 2011 and 2014 as compared to the 
values from 2012, and less than that in the year 2013.® In 
our analysis we ignore these systematical effects as they are 
thought to be unimportant. 


2.5 Data selection criteria and model 
parametrisation 

We are interested in deriving a localised, mean Tatm and 
iatm as a function of month in clear sky conditions. Using 
an unfiltered data, discussed in Sec. 2.2.1 and Sec. 2.4, would 
lead to biased results because only a fraction of the data is 
obtained at the times when there is no cloud cover. The 
fraction of time with clear sky conditions, depends on the 
location and the time of year. In particular, in the lowlands 
surrounding the RT32 observing site the fraction of days 
with the mean daily cloud cover < 20% ranges from ~ 9% 
to ^ 11%, on average (Wos 2010). 

As discussed in Sec. 1, there is an anti-correlation be¬ 
tween RH and solar irradiance, which could statistically hint 
on clear sky conditions, if one selects data by low, ground- 
level RHs. Observations indicate however, that during the 
cloudiest months such selection criterion returns a consid¬ 
erable false-positive rates, as discussed in Sec. 4.6. How¬ 
ever, there is another complication regarding this selection 
scheme. The correlation between water vapour content at 
the ground level and that at higher altitudes is rather weak, 
or non-existent. Some correlation exists only between the 
neighbouring pressure levels in the lower troposphere. For 
this reason, the parametrisation of the selection criterion by 
low RHs at a single, fixed pressure level lacks a very de¬ 
sired feature: that the mean column PWV should decrease 
monotonically as one straightens the selection criterion. 


° In the end of 2011 the station was renewed, hence in 2012 it 
provided the least biased readouts. Therefore, the previous and 
the following years are biassed to a greater extent. 
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Since we are interested in an assessment of the bright¬ 
ness temperature in the conditions that are statistically com¬ 
patible with a clear sky situation in terms of the column 
PWV content, we will filter the data coherently at all pres¬ 
sure levels. For a given month and year, we use only the 
IGRA records that yield 

RH < (3RH,i,j(pH), (5) 

where is the quantile function of the RH distri¬ 

bution at f’th pressure level in j’th month. We use the same 
pn parameter consistently for every pressure level, month, 
and for the ground-based meteorological data. The result¬ 
ing vertical profile will be, in most cases, a combination of 
data taken at different days, but the parametrisation as¬ 
sures that imposing a stronger selection criterion (i.e. low¬ 
ering ph) results in a lower column PWV, as expected. In 
order to average over the year-to-year variability of PWV, 
and approximate the climate for a given month, we split 
the IGRA records into yearly sub-samples, which we anal¬ 
yse separately. Then, for any given month and pressure level 
the mean RH is calculated for all years. The exact value of 
the ph selection parameter, that would be compatible with 
cloudless skies, and would not concern only the driest con¬ 
ditions, is unknown, but can be constrained by independent 
PWV measurements, performed in clear sky conditions (see 
Sec. 2.6). 

The downside of such parametrisation is that the prob¬ 
ability distribution function (PDF) for the Ph selection pa¬ 
rameter, by construction, is zero outside of [0,1] range. This 
is because pn is associated with the probability of obtain¬ 
ing a RH measurement smaller than the quantile function 
(Eq. 5) for a given pressure level, month and year. The final 
monthly-average profile is a multi-year mean, therefore even 
for the maximum value of ph = 1 - i.e. when no data is 
Hltered out - the resulting average model will not be able 
to describe an individual day with PWV values above the 
mean. However, in the current approach we are focused on 
a parametrisation that is suitable for the average clear sky 
model, and will accept these limitations since we are not 
going to analyse the significance of deviations of individual 
measurements from the mean. 

In order to account for the large diurnal variation of 
RH and to enable a meaningful comparison between differ¬ 
ent days, we select the RH data records acquired between 
hours 10 and 14 (UTC-l-1), where the temperatures should 
be least affected by the diurnal variation of the solar irradi- 
ance. Hence, the differences in RH between different days, 
better reflect the actual PWV content variations, and not 
temperature variations. 

2.6 PWV measurements 

We have been monitoring PWV, expressed in column mil¬ 
limetres of water (w), in the clear sky conditions since June 
2013. A compilation from the first data release is sum¬ 
marised in Table 1. The measurements were performed using 
a hand-held MICROTOPS II sun photometer at the Meteo¬ 
rological Observatory of the Department of Meteorology and 
Climatology of the Nicolaus Copernicus University in Torun 
(hereafter DMC). In what follows, we will refer to our sun 
photometer PWV measurements as TR-PWV-DR-1. 

The atmospheric column PWV is calculated along LOS 


towards the Sun, based on photometric measurements of 
water absorption peak at 936 nm, and of the continuum at 
1020 nm (without water absorption). The PWV resulting 
from Beer-Lambert-Bouguer’s law is converted to the zenith 
column water abundance based on the Sun’s zenith distance 
at the time and location of the measurement. The accuracy 
of the measurement is 0.1 mm. Using a reference MICRO¬ 
TOPS H sun photometer we verified that with the factory 
settings the systematical differences to measure PWV be¬ 
tween different instruments is Aw <0.1 mm. 

The data have been collected only in clear sky condi¬ 
tions, which allows us to quantify the PWV variations due 
to air masses carrying different amounts of water vapour, de¬ 
pending on seasons and winds. The data are typically taken 
every hour and every observation contains several measure¬ 
ments that are averaged. In stable clear sky conditions the 
differences between hourly samples are small (relative to our 
measurement precision) and we use a daily mean as a single 
data record. Then, the monthly average, standard deviation 
and extreme values are calculated from daily records (Ta¬ 
ble 1). 

We also utilise the publicly available PWV data from 
AERONET^ robotic stations, which automatically trace the 
Sun and measure the atmospheric water lines spectrum. We 
use the data from the Belsk station in Poland (near War¬ 
saw), collected between (and including) years 2002 and 2014. 
Since the AERONET data are collected automatically and 
exclude only the periods of rain, the raw data may be con¬ 
taminated by the presence of clouds and therefore should be 
more compatible with the average PWV abundances recon¬ 
structed from climatology data: i.e. without imposing any 
selection criteria. The AERONET data are essentially pro¬ 
vided at three levels of post-processing. Level 1.0 is raw data 
(Fig. 5). Level 1.5 data are automatically screened for clouds 
based on number of simple criteria involving time and spec¬ 
tral stability of the measured optical depths (Smirnov et al. 
2000 ), and level 2.0 data, which we use in our analysis, are 
further corrected manually for glitches and any other abnor¬ 
malities that could elude automatic screening. 

The difference between our measurements and the 
AERONET level 2.0 data (Fig. 5) is that the latter have 
clouds removed directly from the LOS towards the Sun, thus 
not guaranteeing a cloudless sky. Our sample however, in 
most cases was taken during cloudless days, at the penalty 
of fewer measurements. 


3 SOLVING RADIATIVE TRANSFER 
THROUGH ATMOSPHERIG LAYERS 

In this section we briefly review the principal equations of 
radiative transfer and outline the calculation scheme. 

Based on a model of the vertical structure of the at¬ 
mosphere we calculate the column PWV (w), atmospheric 
brightness temperature (Tatm), optical depth (r) and the 
corresponding transmittance (tatm = e“’^) using the AM 
program (version 7.2) - a publicly-available radiative trans¬ 
fer solver, developed at the Harvard-Smithsonian Center for 
Astrophysics (Paine 2012). 

^ http://aeronet.gsfc.nasa.gov 
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Propagation of radiation through a medium is given by 
the radiative transfer equation 

^ = ( 6 ) 

as 

where Ii, is the specific intensity and k^, = — is the total 
opacity due to absorption and scattering, is the optical 
depth, and ti, is the emissivity along the propagation path 
s. In the case of local thermodynamic equilibrium (LTE) the 
radiative transfer equation can be rewritten (Wilson et al. 
2009) as: 


1 d/,/ 
Kiy ds 


dh _ j 
dr 


BCT), 


(7) 


where BviT) is the Planck radiance. The solution is given 
by: 

74s) =740)6-""+ f ^ 'B4T(r))e-"‘'dT. (8) 

Jo 

For a thin atmospheric layer, inside which the temperature 
can be assumed constant, Eq. 8 can be integrated as 

74s) = 740)e-""('’> + S4r)(l - (9) 


where the first term on the right hand side is the incident 
radiance that is exponentially attenuated along the propa¬ 
gation path. The second term amounts to the thermal emis¬ 
sion of the layer corrected for the self-absorption, required 
to maintain the assumed LTE. The spectrum of the opti¬ 
cal depth is medium dependent and can be derived from 
the quantum-mechanical properties of atoms and molecules 
present in the medium. The effects leading to a violation 
of the LTE, such as thermal conduction or convection are 
reasonably neglected. 

The radiative transfer equation is solved for each ab¬ 
sorbing (emitting) species and for each atmospheric layer. 
We use Nl ~ 300 stacked layers, each characterised by its 
pressure (P), temperature (T), geometrical thickness, and 
chemical composition: a mixture of nitrogen, oxygen, ozone, 
and water vapour defined in terms of VMRs. The layers 
are homogeneously distributed in log-altitude space, which 
roughly corresponds to a homogeneous distribution in pres¬ 
sure space. We use the same definition of layers for all atmo¬ 
spheric species. The lowermost layer altitude, 2 = 133 m, is 
chosen to coincide with the altitude of RT32. The uppermost 
layer altitude is assumed to be 2max = 60 km. 

It is assumed that the incident (background) radia¬ 
tion has initially Planckian distribution {Bv{T)) with T = 
Tcmb = 2.726 (Fixsen 2009) - the thermodynamic tempera¬ 
ture of the Cosmic Microwave Background radiation (CMB). 
Thus, our definition of Tatm includes the contribution of 
CMB. 

Calculations are performed towards zenith or at direc¬ 
tions located at the zenith distance Zd, and layers are as¬ 
sumed to have a flat geometry. The brightness temperature 
Tatm is obtained by solving 

P,atm = Bu{Ta,tm) (10) 


where P.atm is the output radiance at the lowermost atmo¬ 
spheric layer. The brightness temperature spectrum is cal¬ 
culated within the frequency range v = \1, 60] GHz with 50 
MHz resolution. Details of the physical processes taken into 
account are described in the AM program technical memo 
(Paine 2012). 


4 RESULTS 

4.1 Precipitable Water Vapour 

A compilation of the TR-PWV-DR-1 measurements 
(Sec. 2.6) is shown in Fig. 5, and the corresponding data 
are gathered in Table 1. In this figure solid (cyan) lines rep¬ 
resent the PWV models calculated using the climatological 
data discussed in Sec. 2.1- 2.4. The line width is increased 
with the increases of the pn selection parameter (Sec. 2.5). 
For any given month, the reconstructed column PWV is as¬ 
sumed to occur at mid-month. 

We constrain the value of the pn parameter, which im¬ 
plies the PWV level that is expected to occur in clear sky 
conditions. We use a minimisation, neglecting the cross¬ 
month covariance: 

12 

= ^[{w)i- Mi{pK)\^ /af ( 11 ) 

4 = 1 


where {'w)i and are the i’th month PWV mean and 
standard deviation respectively, and is the model 

PWV value interpolated at the locus of the i’th month data 
point. When the monthly variance estimate is unknown (due 
to missing data), it is interpolated from the neighbouring 
months (with continuity across Dec/Jan). 

Given the PWV data T>, and the reconstructed PWV 
model M , parametrised by pn, we define the posterior prob¬ 
ability for ph using the Bayes theorem: 


T(phIM.T') oc 7:(T'|M,pH)n(pH|M) (12) 


where £(I1|A4 ,Ph) is the likelihood of the data given the 
model and n(pH|A4) is the prior imposed on the parame¬ 
ter probability distribution function (PDF). By design, the 
PDF for the pn parameter is zero outside [0,1] range (see 
Sec. 2.5). We assume the AERONET data as a prior and ob¬ 
tain the maximum posterior constraint on the pn parameter 
and calculate the 68% confidence interval. The likelihood 
function is probed using the sun photometer data (Fig. 6). 
However, we also analyse each of the datasets alone. Con¬ 
straints on the Ph parameter are gathered in Table 2. As 
indicated by the vertical lines in Fig. 6, TR-PWV-DR-1 
and the AERONET/Belsk datasets are compatible at the 
68% CL, but the maximum-likelihood pn parameter value 
is lower for the earlier data, which we attribute to the fact 
that the TR-PWV-DR-1 sample was collected under excel¬ 
lent weather conditions, which guaranteed a cloudless sky, 
and therefore statistically favoured lower PWV abundances. 

From Fig. 5 it is clear that the year-around distribution 
of AERONET/Belsk PWV follows very closely the shape 
of the reconstructed models with high pH values (corre¬ 
sponding to weak data selection). The strength of the cor¬ 
relation is also reflected in the smallness of the best fit 
X^/DOF value, although the dispersion of the individual 
PWV measurements (from which the monthly variance is 
calculated) is relatively large in both data samples. Clearly, 
AERONET/Belsk level 1.0 dataset corresponds to a greater 
Ph value than level 2.0 dataset, as expected. 
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- 1 ^ AERONET level 2.0 (Belsk) 2002-2014 

Sun photometer (Torun) 2013-2015 
Sun photometer (Torun) daily averages 


6 7 

month 


Figure 5. A compilation of PWV measurements. The monthly averages from our TR-PWV-DR-1 sample are marked with blue squares 
while the daily averages are marked with red triangles. The error bars correspond to itlcr deviation from the mean. The black diamonds 
(dots) are the AERONET data level 2.0 (1.0) from the Belsk station. The cyan lines represent our PWV models reconstructed from 
climatology data described in Sec. 2. The thickness of an individual solid line increases with pn- The increases correspond to a progressively 
less aggressive selection of the climatology data by low RH values (Sec. 2.5). 


4.2 Average clear sky atmospheric brightness 
temperature and optical depth: case for 
Torun, Poland 

With the reconstructed vertical structure of the atmosphere 
(Sec. 2) we now determine the local Tatm and r using the 
radiative transfer solver described in Sec. 3. 

In Fig. 7 we plot the best-fit local Tatm models up to 
the Q-band frequencies. The most prominent atmospheric 
feature is associated with the water line that is Doppler- 
and pressure-broadened about the resonance frequency vo « 
22.235 GHz. On the top of the broad band emission a very 
weak and sharp contribution is found due to a thermally- 
excited stimulated emission. This is seen as a tiny increase of 
ATatm ss 0.15 K exactly at the resonance frequency vo, how¬ 
ever effects of fluctuating PWV levels due to atmospheric 
turbulence (see Sect. 4.4) generate Tatm variations of the 
same order in the time scale of hours or less. 

Clearly, Tatm and r are strongly season-dependent, as is 
the level of PWV. The high-frequency tail of Tatm in Fig. 7 
is caused by the contribution from the oxygen O 2 line that 
is insensitive to PWV content. The only way to mitigate 
this emission is to observe at higher altitudes, i.e. through a 
thinner atmosphere. In Fig. 7 the legend provides the value 
of column PWV in millimetres for each model. For compar¬ 
ison, the mean PWV level on the South Pole in the austral 
winter is about 0.26 mm (Bussmann et al. 2005) which is 
~ 25 times lower as compared to the mean PWV content in 
January in Legionowo (Poland) and ~ 12 times lower than 


the mean PWV content in clear sky conditions measured 
in Torun (Table 1). However, since the WV brightness tem¬ 
perature spectrum depends on temperature, a comparison 
of radiative properties between the two sites based on the 
PWV levels difference is not straightforward. 

The impact of ozone is seen only in spectral lines, which 
are weak in the Q-band, as compared to the effects caused by 
atmospheric instabilities. We re-calculate the spectra with 
the resolution of 5 kHz and observe that the lines amplitude 
is well below 1 K above the continuum at the resonance 
frequencies. However, the ozone contribution becomes more 
important at higher frequencies. 

From Fig. 7 it is clear that lower PWV levels result in 
lower Tatm values. The stratospheric PWV levels, as traced 
by satellites, are roughly constant throughout the year. We 
observe that above 13 km the levels are very low, in the 
order of a few /r-meters: w ~ 0.005 mm in January and 
July, and therefore the stratospheric PWV impact on Tatm 
seems unimportant for the ground level sites. 

We use the best fit model determined by the combined 
AERONET/Belsk and Torun data samples (Table 2) to cal¬ 
culate the local Tatm and r for each month. The result is 
shown in Table 3. In that table the error bars correspond to 
the variation of the pn selection parameter within the 68% 
CR. 
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Ph 


Figure 6. Constraints on pn selection parameter (see Sec. 2.5). 
The vertical lines in the bottom panel encompass the 68% CR, 
as summarised in Table 2. The three PDFs are obtained from the 
AERONET/Belsk PWV data (black, dash-dotted), Toruh MI¬ 
CROTOPS measurements (red, solid thin), and the two datasets 
combined (solid thick). 


4.3 Atmospheric approximations and zenith 
distance dependence 

In common observational practice simplifying approxima¬ 
tions are introduced at the cost of accuracy. These approx¬ 
imations rely on assumptions that the atmosphere: (i) is 
composed of a single homogeneous layer, (ii) has a flat ge¬ 
ometry and (iii) is optically thin. In the following sections 
we briefly discuss these assumptions in the light of the re¬ 
constructed atmospheric model. 


4-3.1 Single-layer atmosphere 

A single layer atmosphere has a single physical temperature 
T (Eq. 9) that is often assumed to be somewhere between 
250 K and 290 K. However, given the reconstructed atmo¬ 
spheric model it is easy to derive this value. 

From Eq. 8 the atmospheric emission of a multi-layer, 
flat atmosphere with a vertical temperature profile can be 
written as: 

N 

/..atm = (13) 

where Ti is the i’th layer temperature, is its optical depth, 


Table 1. Compilation of column PWV measurements (TR- 
PWV-DR-1) performed in clear sky conditions using MICRO¬ 
TOPS II sun photometer (Sec. 2.6). 


Location^ 

Period 

Duration^ 

18° 34' 04.8" E, 53° 01' 12.0" N 
2013/06/07 - 2015/04/13 
~ 17 months 


Month 

{U)>° 

Min/Max'* 

# of days® 

Cmt/ 


[mm] 

[mm] 



1 

3.2 ± 1.6 

1.4/4.6 

3 


2 

5.2 ± 1.7 

2.0/7.7 

10 


3 

6.3 ±1.6 

4.4/7.9 

8 


4 

7.8±0.1 

7.8/7.8 

1 

SD 

5 

12.6 ±0.1 

12.6/12.6 

1 

SD 

6 

14.8 ±6.1 

9.3/23.0 

5 


7 

18.5 ±3.5 

16.0/21.0 

2 


8 



0 

ND 

9 



0 

ND 

10 

8.3 ±4.2 

5.6/15.8 

5 


11 



0 

ND 

12 

3.9 ±0.1 

3.9/3.9 

1 

SD 


“Geodetic coordinates of the observation site (DMC). 

^Effective number of observed months. No-observation periods: 
Aug-Sep 2013 and Jul-Oct 2014. 

“Monthly average and a standard deviation of PWV, both calcu¬ 
lated using daily means. 

‘^PWV data extremal values 

“The total number of days observed in a given month. 

^ Comment: SD - single day mean and standard deviation, ND - 
no data 


Table 2. 68% CL constraints on pjj selection parameter from 
the sun photometer data (Table 1), from the AERONET/Belsk 
data and from the two datasets combined. 



Sun Phot. 

AERONET 

Combined 


(Toruh) 

(Belsk, lev. 2.0) 



~ 0.31 

0.04 

NA 

Ph 

(0.28,0.76) 

(0.55,1.00) 

(0.42,0.84) 

MOD(ph) 

^•"±y_0.2i 

0.87t°;i2 

0 62"*”®'^^ 

EX(ph) 

0.53 

0.65 

0.62 


and N is the total number of layers (Sec. 3). is the cal¬ 

culated atmospheric specific intensity, for which the bright¬ 
ness temperature is given by Eq. 10. and is tabulated along 
with T in Table 3. An obvious consequence of considering 
a multi-layer atmosphere is that it does not have a single 
physical temperature, because it is a combination of multiple 
components, each having a different temperature (Eq. 13). 
A multi-layer atmosphere however can be assigned a single¬ 
layer atmosphere equivalent temperature, Tsl, that yields: 


-bB.(rsL(!/))(l (14) 

ing Eq. 13 and Eq. 14 gives 

"" (15) 


(C = 


D = 


ks In//’ 

2hC 


1 - 
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Table 3. 68 % CL constraints on the clear sky, zenith, mean atmospheric brightness temperature (Tatm) and optical depth (r). The values 
for the selected frequencies and months are calculated according to the best fit PWV model selected by the combined AERONET/Belsk 
and Toruh data samples (Table 2). The asterisk (*) marks the errors that are rounded up to appear non-zero at the provided accuracy. 
A and B are the linear-fit constants of the Tatm = A{w) B and tatm = A{w) H- B scaling relations. A and B are written in scientific 
notation with the decimal exponent in parentheses. 



V [GHz] 

5 

15 

22 

30 

5 

15 

22 

30 

Month 

(w) 


Tatm 



T 




[mm] 

[K] 

[K] 

[K] 

[K] 

xlO”® 

xlO -2 

xlO -2 

xlO -2 

1 


r: 1Q+0.01 
o.lU-o 01 

'^•^-0.1 

14.6l[-® 

11 7+®-^ 
-*--*-•'-0.5 

y.ou_o 03 

1 5 + °-^ 
-*^•^-0.1 

4 7+0.8 
' -0.7 

'^•^-0.2 

2 

4 7+0.9 

' -0.9 

o.uy_o 01 

^•'^—0.1 

14.4+;® 

11 7+°-^ 

' -0.4 

y-^o_0.04 

1 4+0.1* 
-*-•^-0.1* 

^•'^-0.6 

'^•^-0.2 

3 



^•^-0.1 

15.5++ 

11.91°;® 

9 09+0.04 

-*^•^-0.1* 

tv 

3 0+0.2 
'^'^-0.2 

4 

7.8tli 

c: 19+0.01 

o.l^-O.oi 

n 7+0.2 
'^•'-0.1 

19.91+ 

I3.0l°;® 

9 34+0.06 

1 ^+^-1 
-*^•^-0.1 

0 0+0.9 
^•'^-0.8 

^•'^-0.2 

5 

11.8+:® 

K 1C+0.01 

o.lO_o 02 

7 o-)-0.2 
' •'^-0.2 

26.91+ 

14.91°-; 

9 30+0.06 

i^.ou_o 07 

1 7+0.1 
' -0.1 

9.3I+ 

4 0+0.3 
^•^-0.3 

6 

16.3+[;^ 

^•^U_o 02 

7 q+0.2 
' ‘^-0.2 

34.71+ 

17 0+°-® 

-*■' '^-0.9 


-^•^-0.1 

12.31+ 

^•^-0.3 

7 

20.9t[® 

c 9 c:+0.02 

o.^O-o 02 

Q c+0.2 

O 2 

42.81+ 

iy.z_o 8 

-0.07 

^•^-0.1 

15.51+ 

'^•-*^-0.3 

8 

18-+?:8 

r 90+0.02 

o.zo_0 02 

Q 1 +0.3 
^•^-0.2 

38.01+ 

‘ -^-0.8 

9 43+0.08 

^•'^-o.l 

13.61+ 

r 7+0.4 

^•'-0.3 

9 

13.0+[-7 

c 1 Q+0.02 

7 k+0.2 
' ‘^-0.2 

29.0I+ 

15.51°;® 

9 39+0.07 

1 Q + 0.1 

^•°-o.i 

10 . 1 +; 

^•^-0.3 

10 

9.3l[® 

K 1 c+O.Ol 
^•10-0.02 

7 q+0.2 
' ‘^-0.2 

22 . 61 + 

13.81°;® 

9 30+0.07 

i^.ou_o 07 

-*^•’^-0.1 

7 7+1.3 
'•'-1.2 

^•^-0.3 

11 

8.7l['® 

c: 10+O.OI 

n 9+0.2 
’^•^-0.2 

21.31+ 

-,0 rr+O.e 

J^'5.0-0.6 

9 40^^'®® 
^^•^^_0.06 

-‘^•’^-0.1 

7 9+0.9 
' -^-1.0 

^•^-0.2 

12 

r 7+1.2 
' -1.0 

r: 1-1+0.01 

'^•^-0.1 

I6.2I+ 

1 9 9+0.5 
-*-^'^-0.5 

05 

-*^•^-0.1 

c O+0.8 
D.CS-o 7 

Q 7+0.2 
'^•'-0.2 

Mean 

10.5 

5.15 

7.1 

24.7 

14.4 

9.42 

1.7 

8.5 

4.5 

St.Dev. 

5.5 

0.05 

0.7 

9.7 

2.6 

0.05 

0.2 

3.7 

0.9 

A 


9.812(-03) 

1.339(-01) 

1.753(+00) 

4.621(-01) 

-4.544(-06) 

4.046(-04) 

6.676(-03) 

1.567(-03) 

B 


5.047(4-00) 

5.713(-l-00) 

6.185+00) 

9.492(4-00) 

9.475(-03) 

1.249+02) 

1.458+02) 

2.811+02) 



u [GHz] 


Figure 7. Average atmospheric brightness temperature as a func¬ 
tion of frequency. The upper (red) lines correspond to July and 
lower (blue) lines correspond to January in Toruh. Within each 
group the thick solid line traces Tatm corresponding to the model 
obtained without imposing any selection criteria on the clima¬ 
tology data (ph = 1.0). Thin solid lines are calculated accord¬ 
ing to the best fit models found using the combined (Toruh + 
AERONET/Belsk) PWV data samples. The best fit pjj selection 
parameter values are given in the plot legend and in Table 2. The 
dashed and dash-dotted lines enclose the 68 % confidence region, 
derived using the Toruh PWV data sample alone (TR-PWV-DR- 
1). All curves are calculated for the zenith distance = 0° at 
the RT32 altitude {z = 0.133 km). 


For a single layer atmosphere N — 1 and Tsl = 7i, and it 
is frequency independent. In general, however, Tsl depends 
on frequency because in different layers the temperature T 
(Eq. 13) is weighted by coefficients (1 — that de¬ 

pend on frequency differently than (1 — e~'^) (Eq. 9), and 
hence the frequency dependence does not cancel. Tsl is ob¬ 
viously season dependent too, and it can be readily derived 
from Eq. 15 using Eq. 10, Tum and r estimates from Ta¬ 
ble 3. However, in the optically thin limit, Tsl is a sensi¬ 
tive function of Tatm and r. Tsl is also a weak function of 
zenith distance. At zenith, in the optically thin limit, it is 
a combination of all atmospheric layers, but near the hori¬ 
zon the lowermost atmospheric layers begin to dominate all 
other contributions. Within the flat-atmosphere model, Tsl 
will reach the ground-level atmospheric temperature at the 
horizon, but it will still vary from one season to another. 
In the optically thick limit, Tsl = Tatm and it also loses its 
frequency dependence because in this limit, the geometry is 
reduced to a single layer case. 

In Fig. 8 we plot the calculated Tsl temperatures. In¬ 
dividual curves result from the best fit atmospheric models 
calculated for the selected months. The spectra were calcu¬ 
lated using Eq. 15, but due to the strong sensitivity of Tsl 
to r and Tatm the values cannot be constrained to a preci¬ 
sion better than 0(10) K from the current data.® The dents 
coincide with the resonance frequencies of water and ozone 
lines. 

For any given frequency the single-layer atmosphere 
equivalent temperature (Tsl) should be most affected by 


® For example, a variation of r by 10 ® at 30 GHz in July implies 
a ~ 4 K change in Tsl (Table 3). 
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Figure 8. Single-layer atmosphere equivalent temperature result¬ 
ing from the best fit model obtained from the combined PWV 
data samples (Table 2). The thin (solid) lines from the bottom 
to the top correspond to January, March and May, and the thick 
(dashed) lines from the top to the bottom are calculated for July, 
September and November. 


the physical temperatures of layers that are most respon¬ 
sible for absorption of radiation at that frequency. About 
50% of total WV coutribution to Tatm takes place above al¬ 
titudes 1-2 km (where the temperature is at least 5-10 K 
below the ground temperatures). At these altitudes there is 
also the biggest concentration of WV. Similarly, the typical 
height scale for ~ 50% of oxygen absorption is about ~ 5 
km with typical temperatures ~ 30 K lower than ground 
temperatures.® Since at 30 GHz the contributions to Tatm 
from oxygen and WV are similar, the expected single-layer 
atmosphere equivalent temperature should be about ~ 20 K 
lower than the ground temperature. The result of a numeri¬ 
cal calculation at 30 GHz is shown in Fig. 9. Clearly, for the 
location of TCfA a reasonable estimate of Tsl, at 30 GHz, 
can be found from the ground level temperatures, offset by 
~ 15.7 K. The offset for other frequencies can be inferred 
from Fig. 8. Given the sensitivity of Tsl (Eq. 15) to r, the 
offset ground-level atmospheric temperature may provide a 
better constraint than measurements of r and Tatm- 

As mentioned earlier, Tsl also depends on zenith dis¬ 
tance, but we find that the dependence is rather weak. Up 
to Zd = 75°, Tsl (2 d) is roughly constant and fixed at its 
zenith value, with accuracy better than {0.2, 0.3, 1.5, 0.8} K 
at frequencies v = {5,15, 22, 30} GHz in July. Deviations 
from the zenith values are even smaller in January. 


Rough estimates found using the best-fit model in July. 



1 2 3 4 5 6 7 8 9 10 11 12 


Month 


Figure 9. Relationship between the single-layer atmosphere 
equivalent temperature, Tsl, at 30 GHz, at the zenith (see Eq. 15) 
and the ground-level month-median atmospheric temperature at 
TCfA. The difference between the two is consistent, within Tier 
error bars, with a constant offset of about 15.7 K. The error 
bars represent the monthly ground-level temperature dispersion 
at TCfA. The median and standard deviation of the ground-level 
atmospheric temperatures are calculated using the meteorologi¬ 
cal data introduced in Sec. 2.4. The error bars are large as they 
include diurnal temperature variations. 

4-S.2 Planar atmosphere 

Radio telescopes with an Az-El mount cannot efficiently ob¬ 
serve sources towards the zenith, thus from practical reasons 
the dependence of Tatm on zenith distance (zd) is important. 
Based on the Bemporad’s air-mass-zenith-distance fitting 
formula (see Schoenberg (1929) or Wilson et al. (2009)), it 
can be seen that the assumption of flatness of the atmo¬ 
sphere - i.e. the atmosphere being composed of flat layers 
stacked one over another - is consistent with sec(2d) scaling 
to within 0.25% (3%) up to Zd < 60° (80°). Within this ap¬ 
proximation (which we use) the atmospheric optical depth 
is given by 

r{zd) = r(0) sec{zd), (16) 

where t(0) = r{zd = 0°). We assess the accuracy of this 
approximation by means of radiative transfer through at¬ 
mospheric layers, whose thicknesses are increased to match 
ray path lengths travelling through spherical layers at given 
zenith angle. As before, we use the same setup of A = 300, 
layers distributed as discussed in Sec. 3, but in this case each 
layer is assumed to be located between geocentric radii Rt 
and Rj and hence its thickness hij = Rj — Ri at the zenith 
angle Zd (measured from the surface of the Earth), is given 
by: 

h,j{zd) = (R? + R] - 2{Rl P CiCj) siC{zd)y^^{17) 


Ci = ^R2csc2(zd)-R| 

Cj = ^R]csc^{zd)-Rl, 

where we assumed spherical Earth. 

For each month, we use the best-fit atmospheric pro¬ 
file and calculate T(zd) relations for a range of zenith dis¬ 
tances, at the frequencies !z[GHz] = {5,15, 22, 30}. Then, we 
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Figure 10. Mean zenith distance dependence of r/r(0°) 
(black,dotted), and (blue,solid and cyan,dash- 

dotted) for the best-fit local atmospheric model. Error bars repre¬ 
sent ±1(T deviation from the mean due to frequency and monthly 
dependence. The short (long) dashed lines trace the season aver¬ 
aged dependence at 5 GHz (22 GHz). 


five to increases of antenna temperature. The relations are 
derived according to the best fit atmospheric model (last 
column of Table 2). The r/r(0) relation for the planar 
case is given by sec(zd) and is not plotted, whereas the 
average t/t(0) for the non-planar case is given by Eq. 18 
(black/dotted line in Fig. 10). The zenith distance depen¬ 
dence of TLra(^d)/'r^tin(0) is given by the following fitting 
formula, which approximates the month and frequency av¬ 
eraged relation below Zd = 89° (dash-dotted/cyan line in 
Fig. 10): 


nuo) / 


1.022 sec(zd) - 3.203 • 10"^ sec^(zd) (19) 


-h9.625 • 10“'‘ sec^(zd) - 1.059 • 10"® sec'‘( 2 d), 

where is atmospheric brightness temperature calcu¬ 

lated for the case when the CMB is not present as a source 
term in the radiative transfer equation, and Tltm(O) = 
TLU^d = 0°). 

Clearly, the difference of the mean Tltm(-Zd)/T’ltm(0) 
from sec(zd) scaling is ~ 0.1 (0.4) or alternatively: (sec( 2 :d) — 
F(tm/7)(tm(0))/sec(zd) « 4% (9%) at Zd = 60° (75°), where 
the averaging is done over twelve months and the four con¬ 
sidered frequencies (Fig. 10). Looking into individual fre¬ 
quencies, the month averaged sec(zd) — r4m/2l(tm(0) ~ 
{0.06,0.09,0.11,0.10} for Zd = 60° at i/ = {5,15,22,30} 
GHz. The effects due to planarity are more than order of 
magnitude smaller. The difference between the averaged 
t/t(0) and Tltm/T'atm(0) scalings increases as the optical 
depth escapes beyond the optically thin limit. This is dis¬ 
cussed in the next section. 


normalise T(zd) at the zenith and average between different 
frequencies and seasons. The differences between individual 
months (January and July) and frequencies vary between 40 
and 53 at the horizon, and between 30 and 35 at Zd = 89°, 
and between 10.6 and 11 at Zd = 85°. 

The following fourth order in sec(zd) fitting function 
provides a good fit to the reconstructed average relation, 
with the maximal deviation below 0.3% at Zd = 89°: 

^ = 1.001 sec(zd) - 3.739 • 10”'^ sec^(zd) (18) 

—4.666 • 10”"^ sec^ (zd) + 6.0366 • 10~® sec"^(2d). 

We do not use Zd = 90° data point to avoid infinite values, 
but since the effects of refraction are not taken into account, 
the usability of the formula is effectively reduced down to 
Zd < 85°, where the dispersion due to a seasonal and fre¬ 
quency dependence becomes similar to uncertainties due to 
neglected refraction. 

According to the fitting formula, the assumption of the 
flat atmosphere (Eq. 16) is accurate to within 0.16 (1.6)% at 
Zd = 60° (80°), where (r(«d)/r(0)) = 1.997±0.001 (5.670± 
0.025) and where the error bars represent the Icr dispersion 
due to season and frequency dependence. 

To address the problem described at the beginning of 
this section, we calculate the average T^tm(zd)/T^tm(0) rela¬ 
tions for the planar and non-planar atmospheres (Fig. 10). 
This is useful because radio telescopes are directly sensi- 



4-S.S Optically thin atmosphere 

Within the RJ approximation B^{T) = 2v^kBTjc?, and 
Eq. 13 becomes: 

N 

Tatm(zd) = H ^ (20) 

The approximation of optically thin atmosphere is often ex¬ 
ploited by utilising the first order expansion in r 

T^Li^d) = rCMB(l -T(Zd)) +TshT{Zd), (21) 

where, as before, we utilised the single layer equivalent tem¬ 
perature at the zenith (see below). We checked that in order 
to visualise the shortcomings of using the optically thin at¬ 
mosphere approximation the differences due to introducing 
the RJ approximation are unimportant. 

Depending on what is being measured, Tcmb contri¬ 
bution is sometimes treated independently from T^tm, as 
if the CMB was not processed by the altitude dependent 
absorption. Below, we quantify this approximation, but in 
general, since the CMB propagates through an absorbing 
and emitting medium, the two contributions (one from the 
atmosphere and the other from the CMB) cannot be treated 
separately. For this reason, our constraints on Tatm include 
the CMB as a source term in Eq. 7. Neglecting the CMB as 
a source term, Eq. 21 for the case of flat atmosphere (Eq. 16) 
becomes: 

riti(^d) = rsLT(0)sec(2d), (22) 
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Figure 11. (Top panel) Accuracy of the atmospheric bright¬ 
ness temperature linear order expansion (Eq. 24) in January 
(blue/thick) and July (red/thin). The best-fit model {pn = 0.62) 
from the combined PWV data samples (Table 2) is used. For 
each frequency the zenith Tsl values are provided in the plot 
legend. {Bottom panel) Accuracy of Tatm ~ Tatm' + Tcmb ap¬ 
proximation expressed as a difference between the actual Tatm 
and the sum of constituents considered separately: Tcmb and at¬ 
mospheric brightness temperature derived without the CMB as a 
source term (T^tm)- 

and what follows is that for a single layer atmosphere 

Tltm(60°) = 2r't^(0°), (23) 

where Ti^tmizd) is the atmospheric brightness temperature 
for Tcmb = 0 K. This approximation is often used in tipping 
scan measurements. 

Increasing the air mass, by looking at greater Zd angles, 
the linear-order approximation in Eq. 22 must fail at high 
values of r as Ts,tm (or T^tm) cannot reach infinity. The ratio 
of Eq. 20 and Eq. 21: 

/./(Tcmb, Tsl, To, Zd) = Tatm/T^tm (24) 

gives a factor by which Tatm is underestimated by using the 
linear-order single-layer approximation. Notice that Eq. 21, 
and consequently Eq. 22, introduce yet another, commonly 
used approximation, that Tsl does not depend on Zd, but 
instead is fixed at its zenith value. As discussed in Sec. 4.3.1, 


this is quite a reasonable approximation for a range of ZdS. 
Tatm in Eq. 24 is calculated numerically within the fiat at¬ 
mosphere approximation, and is calculated according 
to Eq. 21 with T{zd) derived numerically and with Tsl cal¬ 
culated from Eq. 15. 

The factor fi, depends on r, frequency, and vertical 
temperature profile of the atmosphere. Assuming a flat and 
single-layer atmosphere without the CMB, it is easy to see 
that the factor becomes independent of atmospheric tem¬ 
perature and at Zd = 60° amounts to /(Tcmb = 0K,to) = 
(1 _ e-"o=°°("<i))/(rosec(zd)) « 0.965(0.94) at 30 GHz for 
the best-fit value of tq = 0.036(0.061) in January (July) 
(Table 3). 

In the bottom panel of Fig. 11 errors of considering 
Tcmb independently from Tatm are quantified. These are 
the lowest at low frequencies, small optical depths and small 
zenith angles (see plot description for details). 

4-3.4 Relation to observations 

The approximations discussed in the previous sections may 
affect the radio source flux density observations at the levels 
of few to several per-cent. 

In particular, estimates of the atmospheric absorption 
corrections may be biased depending on the assumed ap¬ 
proximations. Flux-density absorption corrections within 
the flat and optically thin atmosphere model are given by: 

St = Sme'^° sec(zrf) ^ sQc{zd)) (25) 

where St is the true flux density and Sm is the mea¬ 
sured flux density. These corrections require an estimate of 
T = T{u,Zd,m), which is a function of frequency {v), zenith 
distance (zd), and month (nr). 

With a single-dish radiometers r is typically estimated 
by measuring the system temperature components at the 
zenith and at Zd = 60° with an implicit assumption of the 
validity of Eq. 23. In the simplest case, when detector linear¬ 
ity is assumed and the spillover, side-lobe and ground pick¬ 
up contributions are neglected, in the RJ approximation, 
the measurement can be defined by the following system of 
linear equations: 

c/E(0°) = Tit,„(0°)-fT,x-fTcMB + AT40°) 

c/E(60°) = A^(60°)Tit^(0°)-fT,x-fTcMB + AT,(60°) 

ef Vibs ~ Tabs Trx, 

(26) 

where Trx is the receiver noise temperature. Tabs is the ab¬ 
sorber temperature and c/ is a conversion factor from the 
voltage, measured at the square-law detector, to the units 
of antenna temperature, and V{zd) is the measured voltage 
at Zd- ATi, is a correction factor weakly dependent on the 
zenith distance (Fig. 11 bottom panel). Usually, at Zd = 60° 
the Am = 2 (and AT^, = 0) is assumed, in accordance with 
the flat atmosphere expectation. 

Substituting thus derived value of T(tm into Eq. 22 with 
an assumed value of Tsl gives an estimate of r. This is a 
possible source of systematical effects. First, as discussed 
in Sec. 4.3.1, Tsl follows ground level temperatures and 
fixing its value will lead to season dependent systematical 
effects. Secondly, Fig. 11 indicates that using an approxi¬ 
mation given by Eq. 22 leads to a slight, but systematical 
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PWV [mm] 


Figure 12. {top panel) Tatm-PWV scaling relations given by 
Eq. 27 and Table 4 (lines). The four sets of points over each 
of the lines represent the best-fit model brightness temperatures 
from Table 3. {bottom panel) Residues between the calculated 
scaling relations described in Sec. 4.4, and their respective fitting 
functions. 


underestimation of r and St- In practice, absorption correc¬ 
tions are applied to both, the target source at and the 
chosen flux calibrator at the zenith distance hence the 
effects of the approximation in Eq. 22 on the relative flux 
density will cancel, as long as the calibrator and the tar¬ 
get source happen to be observed at the same Zd- This is 
however almost never the case, and so the effects of the ap¬ 
proximation will propagate onto the flux density estimates 
through factors ~ (1 -b tq sec(2^''‘^))/(l + tq sec(2:“*)). Since 
the calibrator and target sources are observed at different 
elevations, this effectively leads to an increased variance of 
flux density estimates for individual sources, disparities be¬ 
tween sources of the same intrinsic flux density, and possibly 
to systematical effects depending on source declination with 
respect to the calibrator. 


4.4 Tatm/T-PWV scaling relation 

In clear sky conditions and for a given vertical profile of 
PWV, there is a 1:1 relation between the measured PWV 
and Tatm- In this section we derive Tatm-PWV and r-PWV 
scaling relations for January and July, using the best fit 
model vertical profile estimates from the combined data 
samples (Table 2) . In each profile the PWV content is scaled 


at all altitudes by constant factors, and the corresponding 
Tatm and r are recorded. We calculate the scaling relations 
at frequencies 5,15,22 and 30 GHz and find that the follow¬ 
ing third order polynomial provides a good fit with negligible 
residues: 

3 / X i 




Xi 


w 


1 mm 


(27) 


where X = {Tatm, t} and xt = {at, bt}. The at and bi coeffi¬ 
cients are summarised in Table 4. Terms beyond the linear 
order describe the subtle effects of pressure line broadening 
and self-induced water vapour continuum absorption, both 
non-linearly dependent on H2O VMR (Paine 2015). For r 
the scaling is quadratic and so we use the second order poly¬ 
nomial. In Fig. 12 we visualise Tatm scaling relations and 
the residual errors between these relations and their fitting 
formulas (Eq. 27). 

The fitting formulas defined in Eq. 27 were derived us¬ 
ing the best-fit models for January and July by scaling the 
WV content, hence the Tatm (or r) predicted by these for¬ 
mulas will not exactly match the best-fit model values from 
Table 3. This is because the vertical structure of the atmo¬ 
sphere is month dependent. In order to visualise the ampli¬ 
tude of these differences, for each frequency, we over-plot 
the twelve Tatm values form Table 3 atop the fitting rela¬ 
tion obtained for July (Fig. 12). Tum from different months 
are consistent within Icr error bars quoted in Table 3, but 
a better consistency is reached at high PWV levels, typical 
for summer months, as expected. The year-average PWV- 
Tatm scalings, which result from fitting the data points from 
Fig. 12, are given by A and B coefficients in Table 3. 

Tatm-PWV scaling relation can be used to estimate Tatm 
instabilities due to time varying column PWV. The ampli¬ 
tude of temporal fluctuations of the column PWV, under the 
frozen turbulence hypothesis, should constrain the spatial 
power spectrum of water vapour distribution at the scales 
L, roughly corresponding to VwXt, where At is the time 
interval between subsequent PWV measurements, and v-u, 
is the wind speed at relevant altitudes. Alternatively, on 
the same assumptions, two instantaneous measurements of 
PWV, performed at two locations separated by distance L, 
probe the turbulent spectrum of PWV at the scales roughly 
corresponding to the spatial separation of the observation 
sites. While the spectrum of the atmospheric PWV fluctua¬ 
tions is beyond the scope of this paper, we discuss a rough 
estimate of the amplitude of the Tatm variations in clear sky 
conditions, based on the sun photometer measurements of 
the column PWV. 

Since the Tatm-PWV scaling relation is dominated by 
the linear terms (ao and ai) the ai term can be thought of 
as an approximation of the derivative which quantifies 

Tatm response to variations in PWV (w). Table 4 shows that 
PWV variation of ~ 0.2 mm (~ 0.3 mm) would cause Tatm 
variation of the order ATatm ~ 0.09 (0.14) K at 30 GHz in 
January (July). These are the amplitudes of the PWV varia¬ 
tions in clear sky conditions we actually observe at ~ 1-hour 
time scales. Glearly, low column PWV values also result in 


According to our notation Tatm includes the CMB as a source 
term (Sec. 3). 

The amplitude of the PWV variations quoted are examples 
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Table 4. Tatm-PWV and r-PWV fitting formula coefficients (Eq. 27). The numbers are written in scientific notation with the dec¬ 
imal exponent in parentheses. The ’error’ column contains standard deviations of the difference between the derived scaling and the 
corresponding fitting formula. 


January 


[GHz] 

ao [K] 

ai [K/mm] 

02 [K/mm^l 

03 [K/mm®] 

error [mK] 

bo 

fel [mm ^] 

62 [mm ®] 

5 

5.055(±00) 

8.104(-03) 

7.997(-05) 

-2.841(-09) 

2.441(-06) 

9.350(-03) 

3.108(-05) 

3.047(-07) 

15 

5.738(±00) 

1.282(-01) 

6.849(-04) 

-3.392(-07) 

2.551(-04) 

1.209(-02) 

4.921(-04) 

2.741(-06) 

22 

6.833(±00) 

1.644(4-00) 

-3.827(-03) 

8.033(-06) 

1.434(-01) 

1.650(-02) 

6.391(-03) 

5.896(-06) 

30 

9.477(±00) 

4.603(-01) 

2.396(-03) 

-4.450(-06) 

July 

-2.070(-09) 

6.191(-03) 

2.729(-02) 

1.792(-03) 

1.096(-05) 

5 

5.062(-|-00) 

7.640(-03) 

5.368(-05) 

3.662(-04) 

8.820(-03) 

2.728(-05) 

1.901(-07) 

15 

5.705(±00) 

1.245(-01) 

4.551(-04) 

-2.370(-07) 

3.417(-02) 

1.125(-02) 

4.450(-04) 

1.710(-06) 

22 

6.739(±00) 

1.825(±00) 

-4.879(-03) 

6.163(-06) 

9.216(-01) 

1.515(-02) 

6.616(-03) 

3.679(-06) 

30 

9.230(±00) 

4.455(-01) 

1.555(-03) 

-3.130(-06) 

3.457(-01) 

2.472(-02) 

1.612(-03) 

6.840(-06) 


a more stable atmospheric brightness temperature. By the 
Kolmogorov’s atmospheric turbulence model, the brightness 
fluctuations are characterised by a steep spectrum over a 
wide range of spatial frequencies, and the overall variance is 
dominated by the largest scales. For the wind speed of the 
order 0(1) m/s typical for a calm sky without frontal activ¬ 
ities, 1-hour time scale corresponds to L ~ 0(10) km length 
scales, which coincide with the largest scales from the iner¬ 
tial sub-range, in which the atmospheric turbulence driven 
cascade of kinetic energy transport takes place. Incidentally, 
this rough estimate remains in agreement with theoretical 
predictions from the atmospheric turbulence model (Baars 
2007), which for the same frequency at the sea level gives 
ATatm « 0.0055 Tatm, that is ATatm « {0.06, 0.11} K at 30 
GHz for January and July respectively (Table 3). This is 
also qualitatively compatible with an independent analysis 
involving the near-ground RH variability and wind speed 
measurements (Lew 2016, in preparation). 

4.5 Accuracy limits 

The precision to constrain the local, mean Tatm and r range 
from ~ 0.2% at 5 GHz to ~ 13% at 22 GHz (see Table 3 
for January). These uncertainties refer to the monthly mean 
estimates, rather than to individual PWV measurements, al¬ 
though the dispersion of the individual measurements affects 
the value when model fitting. A random PWV measure¬ 
ment may deviate significantly from the mean even in clear 
sky conditions (e.g. see the scatter from individual PWV 
measurements in Fig. 5 or Table 1). In order to assess the ac¬ 
curacy with which the actual Tum and r can be constrained 
from the statistical analysis of the climatology data, we cal¬ 
culate the impact of the greatest variations, observed in our 
PWV measurements, on the values of Tatm and r. 

The greatest observed variation of PWV in our sample 
occurs in June (Table 1) and ranges from re ~ 9 mm to w ~ 
23 mm and the Icr variation is approximately ±6.1 mm. As¬ 
suming linearity of the PWV-Tatm relation (ai coefficients 
for July in Table 4), this translates to la temperature vari¬ 
ations ATatm < ±{0.05, 0.8,11, 2.7} K at !/ = {5,15, 22, 30} 


from a single day observations in January and July, but the am¬ 
plitude depends on the considered time-scale, and also somewhat 
varies from one day to another. 


GHz, or by about ATatm/Tatm ~ ±{1,9, 26,14} % (see Ta¬ 
ble 3 for July). Clearly, low frequencies are the least sensitive 
to WV variations, as expected. 

4.6 Clear sky detection 

By the analysis of the ground-level solar irradiance (Eq), 
temperature (T) and all-sky images, we observe that in clear 
sky conditions there is a good positive correlation between 
the solar irradiance and air temperature. At such times Eq 
and T are smooth functions of time. On the other hand, 
clouds passing through the LOS towards the Sun generate 
a high frequency noise (Fig. 13 second row panels). We an¬ 
ticipate that this feature can be used for automatic cloud 
detection in meteorological data analyses, but we defer de¬ 
tails to a separate study. 

In Sec. 1 we mentioned that for a fixed pressure level 
the RH has a diurnal variation corresponding to tempera¬ 
ture variations of the atmospheric layer. The RH variation 
anti-correlates with the near-ground atmospheric tempera¬ 
ture and also with the solar irradiance (Fig. 13) detected 
at DMC. However, IGRA data from Legionowo (Sec. 2.2.1) 
are recorded only twice a day (midday and midnight) and 
it is impossible to analyse the stability of diurnal varia¬ 
tions of RH in an attempt to infer the cloud cover. It is 
therefore, interesting to see how feasible it is to select sub¬ 
samples of IGRA data, that would statistically correspond 
to clear sky conditions, given only few observational param¬ 
eters: T,P and Tdew We will investigate this using local sky 
images archived at DMC and T,P and RH measurements 
from TCfA, coarse averaged at 1-hour time scale. The dis¬ 
tance between the two sites is ~ 8 km. 

The amplitude of diurnal temperature variations is the 
greatest in clear sky conditions. This is because the radiative 
cooling of the surface of the Earth is more efficient without 
heat-trapping clouds. Since the RH follows these variabili¬ 
ties in anti-phase, it should be expected that at fixed tem¬ 
peratures a selection based on the lowest daytime RH values 
should statistically correspond to clear sky conditions. On 
the other hand, a thick cloud cover tends to mitigate the 
amplitude of day to night variations in temperature and 
humidity. In this section we report results of a statistical 
analysis aimed to verify the accuracy of this hypothesis. 

The hourly averaged ground-level meteorological data 
(Sec. 2.4) are screened by RH to form a sub-sample of the 
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Table 5. Summary of the sensitivity of the cloud detection al¬ 
gorithm by low ground-level RH values for the two cloud cover 
categories: ’All’ where all types of clouds are considered, and “Low 
and Medium” for which high clouds are treated as a clear sky. 


Cloud cover 



True-Positive fraction [%] 


sky fraction 

Clouds": 


All 


Low & Med. 

(CI“) 

Season^: 

C&H C 

H 

C&H 

C 

H 

< 0.125(1) 


24 

23 

24 

57 

56 

58 

< 0.250(2) 


34 

29 

40 

65 

56 

74 

< 0.375(3) 


48 

35 

62 

76 

62 

89 


“Type of clouds considered. ’Low & Med.’ means that all high 
clouds (if present) were treated as a clear sky. 

*”H’ - hot season: months from April to September, ’C’ - cold 
season: months from October to March, ’C&H’ - all year. 

“Cl - cloud cover index ranging from 0 to 8. For example. Cl = 1 
corresponds to cloud cover sky fraction of 1/8 = 0.125. 

5% driest conditions (sec. 2.4). In order to mitigate effects 
of diurnal temperature variations we consider only the sam¬ 
ples obtained between hours 10:00 and 14:00 of the UTC-l-1 
time. Such a choice roughly corresponds to the times when 
the temperatures should be the most stable (Fig. 13). The 
selection typically picks out two days per each month, de¬ 
pending on data completeness. Next, we visually analyse 
the all-sky images for the selected days and times and as¬ 
sign a mean cloud cover index for two cases: (i) disregarding 
the distinction between high, medium and low clouds, and 
(ii) ignoring high clouds i.e. treating them as clear sky. The 
cloud cover index (Cl) can range from 0 for no clouds sit¬ 
uation, up to 8 for the full sky cloud cover. We disregard 
whether the cloud cover is thin or thick, which is a conser¬ 
vative choice. The result is summarised in Table 5 

Clearly, the selection by RH at the ground level has a 
high false-positive (low true-positive) fraction when all types 
of clouds are considered (see “All” column in Table 5). How¬ 
ever the situation significantly improves when high clouds 
are treated as a clear sky. This result is easy to interpret. 
There is generally a poor correlation between RHs (or the 
corresponding H2O VMRs) at different pressure levels (al¬ 
titudes) hence, there is no guarantee that selection by low 
near-ground RHs will pick out days with no or little high 
cloud cover. Even the straightest selection by RH (5% dri¬ 
est days per month) does not always select fully cloudless 
days. This typically happens for the months during which 
the sky is cloudy most of the time (such as November in 
Poland). An example of a false-positive detection is shown 
in Fig. 13. The two days were selected by the algorithm as 
cloudless. The RH variabilities (bottom panels) and their 
lowest values are similar in both cases. Yet, a visual inspec¬ 
tion revealed that the earlier day was contaminated by high 
clouds and classified as Cl = 6 on the average, as opposed 
to the latter day classified with Cl = 0. It seems however, 
that solar irradiance measured directly to the Sun could be 
effectively used to detect clouds by analysing the high fre¬ 
quency Fourier modes of the time domain signal. When 

We measure the solar irradiance using the CMP22 pyranome- 
ter operating in the spectral range from 200 nm to 3600 nm and 
within ~ 180° FOV. 


3 



Figure 13. Comparison of atmospheric parameters on two differ¬ 
ent days identified as having clear sky conditions, selected on the 
basis of low RH values. The left panels show the correct selection 
(with an average 4-hour cloud cover Cl = 0), and the right pan¬ 
els show an example of a false-positive selection (with an average 
4-hour cloud cover Cl ~ 6). From the top to the bottom the pan¬ 
els show the all-sky camera picture around noon (UTC-l-1), the 
ground-level all-sky solar irradiance, temperature and RH. The 
dashed lines trace the CMP 22 pyranometer sensor temperature. 
While the false-positive selection occurred most likely due to the 
little impact of high clouds on the near-ground RH, the impact 
of such a thin layer of high clouds is clearly significant on the 
detected irradiance. 

treating high clouds as clear sky, however, the selection of 
days with cloud cover below 0.375 is accurate in ~ 76% of 
cases year-round and in ~ 89% of cases during the hot sea¬ 
son (see Table 5). We have not investigated the specificity 
of the estimator - i.e. we have not estimated the fraction of 
clear sky days that eluded detection, but we are aware of 
such cases. 

On average, disregarding the cloud type, the fraction of 
days with cloud cover < 20% is ~ 10% year-round, 12% 
in the hot season^^ and ~ 8% in the cold season for the 
RT32 region (Wos 2010). Based on the data from Table 5 
(interpolated for the cloud cover < 20%), in the case of all 

See Table 5 for the definitions of the hot and cold seasons. 
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types of clouds, using a binomial distribution with the num¬ 
ber of trials fixed by the number of clear sky day candidates 
selected by the algorithm, we estimate that the probabil¬ 
ity of obtaining the reported true-positive rates by chance 
is very low, < 10“'^. Here, we assumed a true-positive out¬ 
come as a successful trial, with the probability defined by 
the aforementioned cloud cover statistics established by the 
local climate. 

While the algorithm of selecting clear sky conditions 
by low, ground-level RHs yields a statistically significant 
result, we have not investigated the true-positive fractions 
where only a single measurement per day is available (case 
for IGRA data). It is possible however, that the algorithm 
can be improved by modifying the time range for data se¬ 
lection (here set between 10:00 and 14:00 of UTC-l-1), or 
by altering the RH threshold (here set at 5% of driest con¬ 
ditions), or by including the long-term pressure variations. 
Due to the reasons described in Sec. 2.5, for the main analy¬ 
sis, we employed external PWV measurements to match the 
data selection criteria to the clear sky requirements. 


5 DISCUSSION 

We used PWV observations (Sec. 2.6) to match the mean 
PWV profile from the sounding data (Sec. 2.2) to the expec¬ 
tations of clear sky conditions. With the current parametri- 
sation (Sec. 2.5), however, we are unable to relate the pn 
best fit value to some other observational quantity that is 
known from meteorological statistics, hence the constrained 
value of the pn selection parameter is potentially useful only 
for the locations with similar climate. Where local clear sky 
PWV measurements are not available, a reasonable approx¬ 
imation of the clear-sky atmospheric WV profile could be 
found using the nearest AERONET data (Sec. 4.1, Fig. 5 
and Fig. 6). 

Although filtering meteorological data by the lowest 5% 
of the ground-level RH, statistically, tends to pick up cloud¬ 
less days (Sec. 4.6 and Table 5), such a single level selection 
does not correspond to the 0.05 value of the pn selection pa¬ 
rameter since, by definition, pn filters the data coherently at 
all altitudes (Sec. 2.5), therefore such an association would 
strongly bias the statistics towards the driest conditions. 


6 CONCLUSIONS 

We use climatological data to reconstruct the vertical struc¬ 
ture of the atmosphere, constrain month-dependent profiles 
of precipitable water vapour (PWV), and predict the atmo¬ 
spheric brightness temperature (Tatm) and optical depth (r) 
at cm-wavelengths. We demonstrate that the nearly global 
coverage of the publicly available climatological data enables 
investigation of almost every location world wide, with the 
spatial resolution of a few hundred kilometres, on average 
(Sec. 2). 

We compare the month-dependence of the column 
PWV, reconstructed for the location of the 32-metre ra¬ 
dio telescope (RT32) located near Toruh (Poland), with the 
AERONET data, collected at the closest station located in 
Belsk. We find that the two are closely correlated through¬ 
out the year, which supports the reliability of the PWV re¬ 


construction from sounding and ground-base meteorological 
data. 

Bearing in mind the prerequisites of radio source con¬ 
tinuum flux density measurements at cm-wavelengths, we 
focus on radiative properties of the atmosphere in clear sky 
conditions. We present a compilation of ~ 17 months of lo¬ 
cal PWV observations collected in clear sky conditions using 
the MICROTOPS H sun photometer (Sec. 4.1). We use these 
observations to devise a selection criterion, which when ap¬ 
plied to the climatological data, enables us to reconstruct 
the vertical structure of the atmosphere that is compatible 
with a cloudless sky. 

Using the reconstructed clear sky PWV profile, and so¬ 
lutions of the radiative transfer through the atmosphere, we 
constrain Tatm and r for the first time for the location of 
RT32 (Sec. 4.2). We also establish PWV-Tatm and PWV-r 
scaling relations (Sec. 4.4) that can be used to constrain at¬ 
mospheric brightness temperature and optical depth in clear 
sky conditions, given an independent estimate of PWV (e.g. 
from a local GPS station). We estimate that in clear sky 
conditions, the mean monthly values of Tatm and r, inferred 
from climatology data, constrain the actual values to within 
±{1,9,26,14} % (at Icr CL) td. u = {5,15,22,30} GHz. 
These estimates should also apply to other locations at sim¬ 
ilar latitudes and a compatible climate (Sec. 4.5). 

We calculate the zenith distance [zd) dependence of 
Tatm and discuss < 10% effects regarding radio-source con¬ 
tinuum flux density measurement calibrations. We discuss 
the implications of using optically thin, single-layer, and 
flat atmosphere approximations in determining the optical 
depth and estimating corrections for atmospheric absorp¬ 
tion (Sec. 4.3.4). For the selected frequencies, we quantify 
deviations of Tatm(2d) and r^Zd) from a simple geometrical 
scaling ~ sec{zd) in the case of non-planar atmosphere. We 
also constrain the physical temperature of the multi-layer 
atmosphere by introducing a single-layer equivalent temper¬ 
ature, which we next connect to the local ground level tem¬ 
peratures by a simple relation (Sec. 4.3). The connection 
should be readily useful when constraining atmospheric op¬ 
tical depth due to absorption and scattering. 

Finally, we discuss the sensitivity of a clear sky selec¬ 
tion criterion involving ground-level relative humidity (RH). 
This criterion can be used to detect cloudless days from data 
that only contain measurements of a few basic atmospheric 
parameters, such as temperature, pressure and dew point, 
which are typically collected by weather balloons and ground 
meteorological stations. By the analysis of archival all-sky 
images, we find that for any given month selecting meteoro¬ 
logical data by the lowest daytime RH can correctly identify 
days with a mean cloud cover below ~ 0.38 in 48% of cases, if 
one disregards whether the sky is obscured by low, medium 
or high level clouds, and in 76% of cases if only low and 
medium level clouds are considered. We find that reproduc¬ 
ing these true-positive fractions by chance is unlikely at odds 
greater than 10"^ : 1 taking into account the local probabil¬ 
ity of cloudless skies, and the number of days used for the 
analysis (Sec. 4.6). We find that the effectiveness of the esti¬ 
mator is increased during the hot season (April-September) 
when the true-positive fraction among the clear-sky days se¬ 
lected by this algorithm reaches 89% (for the case when 
high level clouds are treated as a clear sky). We suspect that 
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the effectiveness of this estimator should be similar in other 
locations with a compatible climate. 
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